# define DISCR_H

# ifndef SCENE_H
# include "scene.h"
# endif

/* ======================================================================================================================== */

/* 2nd order discretization of full "Navier-Stokes Equations" (Momentum Equations) in cylindrical co-ordinate on non-uniform grids. */

/* ---------------Discretizing the "r-momentum" equation---------------------------------------- */
# define RMOMT_TERM_1 ( u[i][j][k]*( (u[i+1][j][k]-u[i][j][k])*(r[i+1]-r[i])/(r[i+2]-r[i+1])\
                               +(u[i][j][k]-u[i-1][j][k])*(r[i+2]-r[i+1])/(r[i+1]-r[i]) )/(r[i+2]-r[i]) )

# define RMOMT_TERM_2 ( (((v[i][j][k]+v[i][j-1][k])*(r[i+2]-r[i+1])+(v[i+1][j][k]+v[i+1][j-1][k])*(r[i+1]-r[i]))/(r[i+1]*(r[i+2]-r[i])*(th[j+2]+th[j+1]-th[j]-th[j-1])))*((u[i][j+1][k]-u[i][j][k])*(th[j+1]-th[j-1])/(th[j+2]-th[j]) + (u[i][j][k]-u[i][j-1][k])*(th[j+2]-th[j])/(th[j+1]-th[j-1])) )

# define RMOMT_TERM_3 ( (((w[i][j][k]+w[i][j][k-1])*(r[i+2]-r[i+1]) + (w[i+1][j][k]+w[i+1][j][k-1])*(r[i+1]-r[i]))/((r[i+2]-r[i])*(z[k+2]+z[k+1]-z[k]-z[k-1])))*((u[i][j][k+1]-u[i][j][k])*(z[k+1]-z[k-1])/(z[k+2]-z[k]) + (u[i][j][k]-u[i][j][k-1])*(z[k+2]-z[k])/(z[k+1]-z[k-1])) )

# define RMOMT_TERM_4 ( pow(0.5*((v[i][j][k]+v[i][j-1][k])*(r[i+2]-r[i+1]) + (v[i+1][j][k]+v[i+1][j-1][k])*(r[i+1]-r[i]))/(r[i+2]-r[i]), 2)/r[i+1] )

#ifndef KVISC_TDEP
# define DEL2U  ( 2*((u[i+1][j][k]-u[i][j][k])/(r[i+2]-r[i+1])-(u[i][j][k]-u[i-1][j][k])/(r[i+1]-r[i]))/(r[i+2]-r[i])\
   +((u[i+1][j][k]-u[i][j][k])*(r[i+1]-r[i])/(r[i+2]-r[i+1])+(u[i][j][k]-u[i-1][j][k])*(r[i+2]-r[i+1])/(r[i+1]-r[i]))/(r[i+1]*(r[i+2]-r[i]))\
   +8*((u[i][j+1][k]-u[i][j][k])/(th[j+2]-th[j])-(u[i][j][k]-u[i][j-1][k])/(th[j+1]-th[j-1]))/(r[i+1]*r[i+1]*(th[j+2]+th[j+1]-th[j]-th[j-1]))\
   +8*((u[i][j][k+1]-u[i][j][k])/(z[k+2]-z[k]) - (u[i][j][k]-u[i][j][k-1])/(z[k+1]-z[k-1]))/(z[k+2]+z[k+1]-z[k]-z[k-1]) )

# define RMOMT_TERM_6  ( (u[i][j][k] + 2*((r[i+2]-r[i+1])*(v[i][j][k]-v[i][j-1][k]) + (r[i+1]-r[i])*(v[i+1][j][k]-v[i+1][j-1][k]))/((r[i+2]-r[i])*(th[j+1]-th[j])))/(r[i+1]*r[i+1]) )

# define RMOMT_VISC  (kvisc0*(DEL2U - RMOMT_TERM_6))

#else

# define RMONT_VISC_1 ( 4*( kvisc(T[i+1][j][k])*(u[i+1][j][k]-u[i][j][k])/(r[i+2]-r[i+1]) - kvisc(T[i][j][k])*(u[i][j][k]-u[i-1][j][k])/(r[i+1]-r[i]) )/(r[i+2]-r[i]) )


# define RMONT_VISC_2 ( ( kvisc(((T[i+1][j][k]*(th[j+2]-th[j+1]) + T[i+1][j+1][k]*(th[j+1]-th[j]))*(r[i+1]-r[i]) + (T[i][j][k]*(th[j+2]-th[j+1])+T[i][j+1][k]*(th[j+1]-th[j]))*(r[i+2]-r[i+1]))/((th[j+2]-th[j])*(r[i+2]-r[i])) )*(2*(v[i+1][j][k]-v[i][j][k])/(r[i+2]-r[i]) - (v[i+1][j][k]*(r[i+1]-r[i]) + v[i][j][k]*(r[i+2]-r[i+1]))/(r[i+1]*(r[i+2]-r[i])) + 2*(u[i][j+1][k]-u[i][j][k])/(r[i+1]*(th[j+2]-th[j]))) - kvisc(((T[i+1][j-1][k]*(th[j+1]-th[j]) + T[i+1][j][k]*(th[j]-th[j-1]))*(r[i+1]-r[i]) + (T[i][j-1][k]*(th[j+1]-th[j])+T[i][j][k]*(th[j]-th[j-1]))*(r[i+2]-r[i+1]))/((th[j+1]-th[j-1])*(r[i+2]-r[i])))*(2*(v[i+1][j-1][k]-v[i][j-1][k])/(r[i+2]-r[i]) - (v[i+1][j-1][k]*(r[i+1]-r[i]) + v[i][j-1][k]*(r[i+2]-r[i+1]))/(r[i+1]*(r[i+2]-r[i])) + 2*(u[i][j][k]-u[i][j-1][k])/(r[i+1]*(th[j+1]-th[j-1]))) )/(r[i+1]*(th[j+1]-th[j])) )


# define RMONT_VISC_3  ( 2*( kvisc(((T[i+1][j][k]*(z[k+2]-z[k+1]) + T[i+1][j][k+1]*(z[k+1]-z[k]))*(r[i+1]-r[i]) + (T[i][j][k]*(z[k+2]-z[k+1])+T[i][j][k+1]*(z[k+1]-z[k]))*(r[i+2]-r[i+1]))/((r[i+2]-r[i])*(z[k+2]-z[k])))*((u[i][j][k+1]-u[i][j][k])/(z[k+2]-z[k])+(w[i+1][j][k]-w[i][j][k])/(r[i+2]-r[i])) - kvisc(((T[i+1][j][k-1]*(z[k+1]-z[k]) + T[i+1][j][k]*(z[k]-z[k-1]))*(r[i+1]-r[i]) + (T[i][j][k-1]*(z[k+1]-z[k])+T[i][j][k]*(z[k]-z[k-1]))*(r[i+2]-r[i+1]))/((r[i+2]-r[i])*(z[k+1]-z[k-1])))*((u[i][j][k]-u[i][j][k-1])/(z[k+1]-z[k-1])+(w[i+1][j][k-1]-w[i][j][k-1])/(r[i+2]-r[i])) )/(z[k+1]-z[k]) )


# define RMONT_VISC_4  ( 2*kvisc( (T[i+1][j][k]*(r[i+1]-r[i])+T[i][j][k]*(r[i+2]-r[i+1]))/(r[i+2]-r[i]) )*( ((r[i+1]-r[i])*(u[i+1][j][k]-u[i][j][k])/(r[i+2]-r[i+1])+(r[i+2]-r[i+1])*(u[i][j][k]-u[i-1][j][k])/(r[i+1]-r[i]))/(r[i+2]-r[i])-(u[i][j][k]+((v[i+1][j][k]-v[i+1][j-1][k])*(r[i+1]-r[i]) + (v[i][j][k]-v[i][j-1][k])*(r[i+2]-r[i+1]))/((th[j+1]-th[j])*(r[i+2]-r[i])))/r[i+1] )/r[i+1] )


# define RMOMT_VISC   (RMONT_VISC_1 + RMONT_VISC_2 + RMONT_VISC_3 + RMONT_VISC_4)
#endif

# define AA (0 - RMOMT_TERM_1 - RMOMT_TERM_2 - RMOMT_TERM_3 + RMOMT_TERM_4 + RMOMT_VISC)
/* --------------------------------------------------------------------------------------- */


/* ---------------Discretizing the "r-momentum" equation at central point---------------------- */
#ifndef KVISC_TDEP

# define rA_CENTER ( kvisc0*(u[i+1][j][k]-u[i-1][j][k])/(r[i+2]-r[i]) )

#else

# define rA_CENTER ( 2*( kvisc(0.5*((T[i][j][k]+T[i+1][j][k])*(th[j+2]-th[j+1]) + (T[i][j+1][k]+T[i+1][j+1][k])*(th[j+1]-th[j]))/(th[j+2]-th[j]))*(v[i+1][j][k]-v[i][j][k]) - kvisc(0.5*((T[i][j-1][k]+T[i+1][j-1][k])*(th[j+1]-th[j]) + (T[i][j][k]+T[i+1][j][k])*(th[j]-th[j-1]))/(th[j+1]-th[j-1]))*(v[i+1][j-1][k]-v[i][j-1][k]) - 2*kvisc(0.5*(T[i+1][j][k]+T[i][j][k]))*(v[i+1][j][k]-v[i][j][k]-v[i+1][j-1][k]+v[i][j-1][k]) )/((r[i+2]-r[i])*(th[j+1]-th[j])) )

#endif

/* --------------------------------------------------------------------------------------------- */


/* -----------------------------Discretizing the "theta-momentum" equation-------------- */
# define THMOMT_TERM_1 ( (((u[i][j][k]+u[i-1][j][k])*(th[j+2]-th[j+1]) + (u[i][j+1][k]+u[i-1][j+1][k])*(th[j+1]-th[j]))/((th[j+2]-th[j])*(r[i+2]+r[i+1]-r[i]-r[i-1])))*((v[i+1][j][k]-v[i][j][k])*(r[i+1]-r[i-1])/(r[i+2]-r[i]) + (v[i][j][k]-v[i-1][j][k])*(r[i+2]-r[i])/(r[i+1]-r[i-1])) )

# define THMOMT_TERM_2 ( 4*v[i][j][k]*((th[j+1]-th[j-1])*(v[i][j+1][k]-v[i][j][k])/(th[j+2]-th[j])+(th[j+2]-th[j])*(v[i][j][k]-v[i][j-1][k])/(th[j+1]-th[j-1]))/((r[i]+r[i+1])*(th[j+2]+th[j+1]-th[j]-th[j-1])) )

# define THMOMT_TERM_3 ( (((w[i][j][k]+w[i][j][k-1])*(th[j+2]-th[j+1]) + (w[i][j+1][k]+w[i][j+1][k-1])*(th[j+1]-th[j]))/((th[j+2]-th[j])*(z[k+2]+z[k+1]-z[k]-z[k-1])))*((v[i][j][k+1]-v[i][j][k])*(z[k+1]-z[k-1])/(z[k+2]-z[k]) + (v[i][j][k]-v[i][j][k-1])*(z[k+2]-z[k])/(z[k+1]-z[k-1])) )

# define THMOMT_TERM_4 ( v[i][j][k]*((u[i][j][k]+u[i-1][j][k])*(th[j+2]-th[j+1])+(u[i][j+1][k]+u[i-1][j+1][k])*(th[j+1]-th[j]))/((th[j+2]-th[j])*(r[i]+r[i+1])) )


#ifndef KVISC_TDEP

# define DEL2V  ( 8*((v[i+1][j][k]-v[i][j][k])/(r[i+2]-r[i])-(v[i][j][k]-v[i-1][j][k])/(r[i+1]-r[i-1]))/(r[i+2]+r[i+1]-r[i]-r[i-1])\
+ 4*((v[i+1][j][k]-v[i][j][k])*(r[i+1]-r[i-1])/(r[i+2]-r[i])+(v[i][j][k]-v[i-1][j][k])*(r[i+2]-r[i])/(r[i+1]-r[i-1]))\
                                                                        /((r[i]+r[i+1])*(r[i+2]+r[i+1]-r[i]-r[i-1]))\
+8*((v[i][j+1][k]-v[i][j][k])/(th[j+2]-th[j+1])-(v[i][j][k]-v[i][j-1][k])/(th[j+1]-th[j]))/((r[i]+r[i+1])*(r[i]+r[i+1])*(th[j+2]-th[j]))\
+ 8*((v[i][j][k+1]-v[i][j][k])/(z[k+2]-z[k])-(v[i][j][k]-v[i][j][k-1])/(z[k+1]-z[k-1]))/(z[k+2]+z[k+1]-z[k]-z[k-1]) )

# define THMOMT_TERM_6 (4*(v[i][j][k] - 2*(u[i][j+1][k]+u[i-1][j+1][k]-u[i][j][k]-u[i-1][j][k])/(th[j+2]-th[j]))/((r[i]+r[i+1])*(r[i]+r[i+1])))

# define  THMOMT_VISC  kvisc0*(DEL2V - THMOMT_TERM_6)

#else

# define THMOMT_VISC_1  ( ( kvisc(((T[i+1][j][k]*(th[j+2]-th[j+1])+T[i+1][j+1][k]*(th[j+1]-th[j]))*(r[i+1]-r[i])+(T[i][j][k]*(th[j+2]-th[j+1])+T[i][j+1][k]*(th[j+1]-th[j]))*(r[i+2]-r[i+1]))/((th[j+2]-th[j])*(r[i+2]-r[i])))*(2*(v[i+1][j][k]-v[i][j][k])/(r[i+2]-r[i])+(-(v[i+1][j][k]*(r[i+1]-r[i])+v[i][j][k]*(r[i+2]-r[i+1]))/(r[i+2]-r[i])+2*(u[i][j+1][k]-u[i][j][k])/(th[j+2]-th[j]))/r[i+1])-kvisc(((T[i][j][k]*(th[j+2]-th[j+1])+T[i][j+1][k]*(th[j+1]-th[j]))*(r[i]-r[i-1])+(T[i-1][j][k]*(th[j+2]-th[j+1])+T[i-1][j+1][k]*(th[j+1]-th[j]))*(r[i+1]-r[i]))/((th[j+2]-th[j])*(r[i+1]-r[i-1])))*(i==NrG?2*(u[i][j+1][k]-u[i-2][j+1][k]-u[i][j][k]+u[i-2][j][k])/((r[i+1]-r[i-1])*(th[j+2]-th[j])):2*(v[i][j][k]-v[i-1][j][k])/(r[i+1]-r[i-1])+(-(v[i][j][k]*(r[i]-r[i-1])+v[i-1][j][k]*(r[i+1]-r[i]))/(r[i+1]-r[i-1])+2*(u[i-1][j+1][k]-u[i-1][j][k])/(th[j+2]-th[j]))/r[i]) )/(r[i+1]-r[i]) )

# define  THMOMT_VISC_2  ( 16*( kvisc(T[i][j+1][k])*((v[i][j+1][k]-v[i][j][k])/(th[j+2]-th[j+1])+0.5*(u[i][j+1][k]+u[i-1][j+1][k]))-kvisc(T[i][j][k])*((v[i][j][k]-v[i][j-1][k])/(th[j+1]-th[j])+0.5*(u[i][j][k]+u[i-1][j][k])) )/((r[i]+r[i+1])*(r[i]+r[i+1])*(th[j+2]-th[j])) )

# define  THMOMT_VISC_3   ( 2*( kvisc(((T[i][j][k]*(th[j+2]-th[j+1])+T[i][j+1][k]*(th[j+1]-th[j]))*(z[k+2]-z[k+1])+(T[i][j][k+1]*(th[j+2]-th[j+1])+T[i][j+1][k+1]*(th[j+1]-th[j]))*(z[k+1]-z[k]))/((th[j+2]-th[j])*(z[k+2]-z[k])))*((v[i][j][k+1]-v[i][j][k])/(z[k+2]-z[k]) + 2*(w[i][j+1][k]-w[i][j][k])/((r[i]+r[i+1])*(th[j+2]-th[j])))-kvisc(((T[i][j][k-1]*(th[j+2]-th[j+1])+T[i][j+1][k-1]*(th[j+1]-th[j]))*(z[k+1]-z[k])+(T[i][j][k]*(th[j+2]-th[j+1])+T[i][j+1][k]*(th[j+1]-th[j]))*(z[k]-z[k-1]))/((th[j+2]-th[j])*(z[k+1]-z[k-1])))*((v[i][j][k]-v[i][j][k-1])/(z[k+1]-z[k-1]) + 2*(w[i][j+1][k-1]-w[i][j][k-1])/((r[i]+r[i+1])*(th[j+2]-th[j]))) )/(z[k+1]-z[k]) )

# define  THMOMT_VISC_4  ( 8*kvisc((T[i][j][k]*(th[j+2]-th[j+1])+T[i][j+1][k]*(th[j+1]-th[j]))/(th[j+2]-th[j]))*(((r[i+1]-r[i-1])*(v[i+1][j][k]-v[i][j][k])/(r[i+2]-r[i])+(r[i+2]-r[i])*(v[i][j][k]-v[i-1][j][k])/(r[i+1]-r[i-1]))/(r[i+2]+r[i+1]-r[i]-r[i-1])-v[i][j][k]/(r[i]+r[i+1]) + (u[i][j+1][k]+u[i-1][j+1][k]-u[i][j][k]-u[i-1][j][k])/((r[i]+r[i+1])*(th[j+2]-th[j])))/(r[i]+r[i+1]) )

# define  THMOMT_VISC   (THMOMT_VISC_1 + THMOMT_VISC_2 + THMOMT_VISC_3 + THMOMT_VISC_4)

#endif

# define BB (0 - THMOMT_TERM_1 - THMOMT_TERM_2 - THMOMT_TERM_3 - THMOMT_TERM_4 + THMOMT_VISC)
/* --------------------------------------------------------------------------------- */


/* -----------------------Discretizing the "z-momentum" equation------------------------- */
# define ZMOMT_TERM_1 ( (((u[i][j][k]+u[i-1][j][k])*(z[k+2]-z[k+1]) + (u[i][j][k+1]+u[i-1][j][k+1])*(z[k+1]-z[k]))/((z[k+2]-z[k])*(r[i+2]+r[i+1]-r[i]-r[i-1])))*((w[i+1][j][k]-w[i][j][k])*(r[i+1]-r[i-1])/(r[i+2]-r[i]) + (w[i][j][k]-w[i-1][j][k])*(r[i+2]-r[i])/(r[i+1]-r[i-1])) )

# define ZMOMT_TERM_2 ( (2*((v[i][j][k]+v[i][j-1][k])*(z[k+2]-z[k+1])+(v[i][j][k+1]+v[i][j-1][k+1])*(z[k+1]-z[k]))/((z[k+2]-z[k])*(r[i]+r[i+1])*(th[j+2]+th[j+1]-th[j]-th[j-1])))*((w[i][j+1][k]-w[i][j][k])*(th[j+1]-th[j-1])/(th[j+2]-th[j])+(w[i][j][k]-w[i][j-1][k])*(th[j+2]-th[j])/(th[j+1]-th[j-1])) )

# define ZMOMT_TERM_3 ( w[i][j][k]*((z[k+1]-z[k])*(w[i][j][k+1]-w[i][j][k])/(z[k+2]-z[k+1])+(z[k+2]-z[k+1])*(w[i][j][k]-w[i][j][k-1])/(z[k+1]-z[k]))/(z[k+2]-z[k]) )


#ifndef KVISC_TDEP

# define DEL2W  ( 8*((w[i+1][j][k]-w[i][j][k])/(r[i+2]-r[i])-(w[i][j][k]-w[i-1][j][k])/(r[i+1]-r[i-1]))/(r[i+2]+r[i+1]-r[i]-r[i-1])\
+4*((w[i+1][j][k]-w[i][j][k])*(r[i+1]-r[i-1])/(r[i+2]-r[i])+(w[i][j][k]-w[i-1][j][k])*(r[i+2]-r[i])/(r[i+1]-r[i-1]))/((r[i]+r[i+1])*(r[i+2]+r[i+1]-r[i]-r[i-1]))\
+32*((w[i][j+1][k]-w[i][j][k])/(th[j+2]-th[j])-(w[i][j][k]-w[i][j-1][k])/(th[j+1]-th[j-1]))/((r[i]+r[i+1])*(r[i]+r[i+1])*(th[j+2]+th[j+1]-th[j]-th[j-1]))\
+2*((w[i][j][k+1]-w[i][j][k])/(z[k+2]-z[k+1])-(w[i][j][k]-w[i][j][k-1])/(z[k+1]-z[k]))/(z[k+2]-z[k]) )

# define ZMOMT_VISC  (kvisc0*DEL2W)

#else

# define ZMOMT_VISC_1   ( 2*( kvisc(((T[i+1][j][k]*(z[k+2]-z[k+1])+T[i+1][j][k+1]*(z[k+1]-z[k]))*(r[i+1]-r[i])+(T[i][j][k]*(z[k+2]-z[k+1])+T[i][j][k+1]*(z[k+1]-z[k]))*(r[i+2]-r[i+1]))/((r[i+2]-r[i])*(z[k+2]-z[k])))*((u[i][j][k+1]-u[i][j][k])/(z[k+2]-z[k])+(w[i+1][j][k]-w[i][j][k])/(r[i+2]-r[i]))-kvisc(((T[i][j][k]*(z[k+2]-z[k+1])+T[i][j][k+1]*(z[k+1]-z[k]))*(r[i]-r[i-1])+(T[i-1][j][k]*(z[k+2]-z[k+1])+T[i-1][j][k+1]*(z[k+1]-z[k]))*(r[i+1]-r[i]))/((r[i+1]-r[i-1])*(z[k+2]-z[k])))*((u[i-1][j][k+1]-u[i-1][j][k])/(z[k+2]-z[k])+(w[i][j][k]-w[i-1][j][k])/(r[i+1]-r[i-1])) )/(r[i+1]-r[i]) )

# define ZMOMT_VISC_2   ( 4*( kvisc(((T[i][j][k]*(th[j+2]-th[j+1])+T[i][j+1][k]*(th[j+1]-th[j]))*(z[k+2]-z[k+1])+(T[i][j][k+1]*(th[j+2]-th[j+1])+T[i][j+1][k+1]*(th[j+1]-th[j]))*(z[k+1]-z[k]))/((th[j+2]-th[j])*(z[k+2]-z[k])))*((v[i][j][k+1]-v[i][j][k])/(z[k+2]-z[k])+2*(w[i][j+1][k]-w[i][j][k])/((r[i]+r[i+1])*(th[j+2]-th[j])))-kvisc(((T[i][j-1][k]*(th[j+1]-th[j])+T[i][j][k]*(th[j]-th[j-1]))*(z[k+2]-z[k+1])+(T[i][j-1][k+1]*(th[j+1]-th[j])+T[i][j][k+1]*(th[j]-th[j-1]))*(z[k+1]-z[k]))/((th[j+1]-th[j-1])*(z[k+2]-z[k])))*((v[i][j-1][k+1]-v[i][j-1][k])/(z[k+2]-z[k])+2*(w[i][j][k]-w[i][j-1][k])/((r[i]+r[i+1])*(th[j+1]-th[j-1]))) )/((th[j+1]-th[j])*(r[i]+r[i+1])) )

# define ZMOMT_VISC_3   ( 4*( kvisc(T[i][j][k+1])*(w[i][j][k+1]-w[i][j][k])/(z[k+2]-z[k+1])-kvisc(T[i][j][k])*(w[i][j][k]-w[i][j][k-1])/(z[k+1]-z[k]) )/(z[k+2]-z[k]) )

# define ZMOMT_VISC_4  ( 2*kvisc((T[i][j][k]*(z[k+2]-z[k+1])+T[i][j][k+1]*(z[k+1]-z[k]))/(z[k+2]-z[k]))*((u[i][j][k+1]+u[i-1][j][k+1]-u[i][j][k]-u[i-1][j][k])/(z[k+2]-z[k])+2*((r[i+1]-r[i-1])*(w[i+1][j][k]-w[i][j][k])/(r[i+2]-r[i])+(r[i+2]-r[i])*(w[i][j][k]-w[i-1][j][k])/(r[i+1]-r[i-1]))/(r[i+2]+r[i+1]-r[i]-r[i-1]) )/(r[i]+r[i+1]) )

# define ZMOMT_VISC  (ZMOMT_VISC_1 + ZMOMT_VISC_2 + ZMOMT_VISC_3 + ZMOMT_VISC_4)

#endif

#ifndef DENS_TDEP
# define ZMOMT_TERM_5  GRAVITY
#else
# define ZMOMT_TERM_5  ( GRAVITY*dens(0.5*(T[i][j][k]+T[i][j][k+1]))/density0 )
#endif

# define CC (0 - ZMOMT_TERM_1 - ZMOMT_TERM_2 - ZMOMT_TERM_3 + ZMOMT_VISC - ZMOMT_TERM_5)
/* ------------------------------------------------------------------------------------------ */


/* ---------------------------- Discretizing for THETA-VORTICITY ------------------------- */
# define VORT_TH_ijk  (\
           0.5*( ((u[i][j][k]+u[i-1][j][k])*(z[k+2]-z[k+1])+(u[i][j][k+1]-u[i-1][j][k+1])*(z[k+1]-z[k]))/(z[k+2]-z[k])\
               - ((u[i][j][k-1]+u[i-1][j][k-1])*(z[k+1]-z[k])+(u[i][j][k]+u[i-1][j][k])*(z[k]-z[k-1]))/(z[k+1]-z[k-1]) ) / (z[k+1]-z[k])\
         - 0.5*( ((w[i][j][k]+w[i][j][k-1])*(r[i+1]-r[i+1])+(w[i+1][j][k]+w[i+1][j][k-1])*(r[i+1]-r[i]))/(r[i+2]-r[i])\
               - ((w[i-1][j][k]+w[i-1][j][k-1])*(r[i+1]-r[i])+(w[i][j][k]+w[i][j][k-1])*(r[i]-r[i-1]))/(r[i+1]-r[i-1]) ) / (r[i+1]-r[i])\
                      )
/* -------------------------------------------------------------------------------------------*/


/* ---------------------------- Discretizing for Divergence (DEL.V) ------------------------- */
# define DIV_ijk (\
          (2/(r[i]+r[i+1]))*( (r[i+1]*u[i][j][k] - r[i]*u[i-1][j][k])/(r[i+1]-r[i]) + (v[i][j][k] - v[i][j-1][k])/(th[j+1]-th[j]) )\
        + (w[i][j][k] - w[i][j][k-1])/(z[k+1]-z[k])\
                 )
/* -------------------------------------------------------------------------------------------*/


/* ========================================================================================================================== */
/* Discretization of div2p used in "vel_advance()". This kind of discretization can be used for any quantity defined on
   CELL-center (e.g. T). However, this CAN'T be used for something defined on FACE-center (e.g. u, v, w). */
# define DIV2P_ijk (\
     8*((p[i+1][j][k] - p[i][j][k])/(r[i+2]-r[i]) - (p[i][j][k]-p[i-1][j][k])/(r[i+1]-r[i-1]))/(r[i+2]+r[i+1]-r[i]-r[i-1])\
   + 4*((p[i+1][j][k] - p[i][j][k])*(r[i+1]-r[i-1])/(r[i+2]-r[i]) + (p[i][j][k] - p[i-1][j][k])*(r[i+2]-r[i])/(r[i+1]-r[i-1]))\
          /((r[i+1]+r[i])*(r[i+2]+r[i+1]-r[i]-r[i-1]))\
   + 32*((p[i][j+1][k] - p[i][j][k])/(th[j+2]-th[j]) - (p[i][j][k] - p[i][j-1][k])/(th[j+1]-th[j-1]))\
          /((r[i+1]+r[i])*(r[i+1]+r[i])*(th[j+2] + th[j+1] - th[j] - th[j-1]))\
   + 8*((p[i][j][k+1] - p[i][j][k])/(z[k+2]-z[k]) - (p[i][j][k] - p[i][j][k-1])/(z[k+1]-z[k-1]))/(z[k+2]+z[k+1]-z[k]-z[k-1])\
                   )
/* ========================================================================================================================== */



/* ======================================================================================================================== */

/* Discretization of "Pressure Poisson Equation". */

# define POISS_BETA_NONUNIF (8*(p[i+1][j][k]/(r[i+2]-r[i]) + p[i-1][j][k]/(r[i+1]-r[i-1]))/(r[i+2]+r[i+1]-r[i]-r[i-1])\
       + (4/(r[i+2]+r[i+1]-r[i]-r[i-1]))*((r[i+1]-r[i-1])*p[i+1][j][k]/(r[i+2]-r[i]) - (r[i+2]-r[i])*p[i-1][j][k]/(r[i+1]-r[i-1]))/(r[i+1]+r[i])\
       + 32*(p[i][j+1][k]/(th[j+2]-th[j]) + p[i][j-1][k]/(th[j+1]-th[j-1]))/((r[i+1]+r[i])*(r[i+1]+r[i])*(th[j+2]+th[j+1]-th[j]-th[j-1]))\
       + 8*(p[i][j][k+1]/(z[k+2]-z[k])+p[i][j][k-1]/(z[k+1]-z[k-1]))/(z[k+2]+z[k+1]-z[k]-z[k-1])\
       - d[i][j][k])
//       - d[i-NrG][j-NthG][k-NzG])

# define POISS_ALPHA_NONUNIF ((8/(r[i+2]+r[i+1]-r[i]-r[i-1]))*(1/(r[i+2]-r[i])+1/(r[i+1]-r[i-1]))\
         - (4/(r[i+2]+r[i+1]-r[i]-r[i-1]))*((r[i+2]-r[i])/(r[i+1]-r[i-1])-(r[i+1]-r[i-1])/(r[i+2]-r[i]))/(r[i+1]+r[i])\
         + 32*(1/(th[j+2]-th[j])+1/(th[j+1]-th[j-1]))/((r[i+1]+r[i])*(r[i+1]+r[i])*(th[j+2]+th[j+1]-th[j]-th[j-1]))\
         + 8*(1/(z[k+2]-z[k])+1/(z[k+1]-z[k-1]))/(z[k+2]+z[k+1]-z[k]-z[k-1]))

# define  P_ijk_NEW  (POISS_BETA_NONUNIF/POISS_ALPHA_NONUNIF)

/* ======================================================================================================================== */


/* Discretization of the Energy Equation. */

#ifdef SOLVE_ENERGY

/* Discretizing the convection term. */
# ifdef EEQ_DISSIP_FACT
#   define dTdr_CD2  ( 2*((r[i+1]-r[i-1])*(T[i+1][j][k]-T[i][j][k])/(r[i+2]-r[i]) + (r[i+2]-r[i])*(T[i][j][k]-T[i-1][j][k])/(r[i+1]-r[i-1]))/(r[i+2]+r[i+1]-r[i]-r[i-1]) )

#   define dTdr_DISSIP ( (T[i+1][j][k]-T[i][j][k])/(r[i+2]-r[i])-(T[i][j][k]-T[i-1][j][k])/(r[i+1]-r[i-1]) )

#   define udTdr_UD1 ( 0.5*(u[i][j][k]+u[i-1][j][k])*(dTdr_CD2\
           + EEQ_DISSIP_FACT*2*((u[i][j][k]+u[i-1][j][k] > 0) ? -r[i+1]+r[i-1]:r[i+2]-r[i])*dTdr_DISSIP/(r[i+2]+r[i+1]-r[i]-r[i-1])) )



#   define dTdth_CD2 ( 2*((th[j+1]-th[j-1])*(T[i][j+1][k]-T[i][j][k])/(th[j+2]-th[j])+(th[j+2]-th[j])*(T[i][j][k]-T[i][j-1][k])/(th[j+1]-th[j-1]))/(th[j+2]+th[j+1]-th[j]-th[j-1]) )

#   define dTdth_DISSIP ( (T[i][j+1][k]-T[i][j][k])/(th[j+2]-th[j])-(T[i][j][k]-T[i][j-1][k])/(th[j+1]-th[j-1]) )

#   define vrdTdth_UD1  ( (v[i][j][k]+v[i][j-1][k])*(dTdth_CD2\
 + EEQ_DISSIP_FACT*2*((v[i][j][k]+v[i][j-1][k] > 0) ? -th[j+1]+th[j-1]:th[j+2]-th[j])*dTdth_DISSIP/(th[j+2]+th[j+1]-th[j]-th[j-1]))/(r[i]+r[i+1]) )



#   define dTdz_CD2  ( 2*((z[k+1]-z[k-1])*(T[i][j][k+1]-T[i][j][k])/(z[k+2]-z[k]) + (z[k+2]-z[k])*(T[i][j][k]-T[i][j][k-1])/(z[k+1]-z[k-1]))/(z[k+2]+z[k+1]-z[k]-z[k-1]) )

#   define dTdz_DISSIP ( (T[i][j][k+1]-T[i][j][k])/(z[k+2]-z[k])-(T[i][j][k]-T[i][j][k-1])/(z[k+1]-z[k-1]) )

#   define wdTdz_UD1 ( 0.5*(w[i][j][k]+w[i][j][k-1])*(dTdz_CD2\
   + EEQ_DISSIP_FACT*2*((w[i][j][k]+w[i][j][k-1] > 0) ? -z[k+1]+z[k-1]:z[k+2]-z[k])*dTdz_DISSIP/(z[k+2]+z[k+1]-z[k]-z[k-1])) )

# else

#   define udTdr_UD1 ( (u[i][j][k]+u[i-1][j][k])*( (u[i][j][k]+u[i-1][j][k] > 0) ? (T[i][j][k]-T[i-1][j][k])/(r[i+1]-r[i-1]) : (T[i+1][j][k]-T[i][j][k])/(r[i+2]-r[i]) ) )

#   define vrdTdth_UD1 ( 2*(v[i][j][k]+v[i][j-1][k])*( (v[i][j][k]+v[i][j-1][k] > 0) ? (T[i][j][k]-T[i][j-1][k])/(th[j+1]-th[j-1]) : (T[i][j+1][k]-T[i][j][k])/(th[j+2]-th[j]) )/(r[i]+r[i+1]) )

#   define wdTdz_UD1  ( (w[i][j][k]+w[i][j][k-1])*( (w[i][j][k] + w[i][j][k-1] > 0) ? (T[i][j][k]-T[i][j][k-1])/(z[k+1]-z[k-1]) : (T[i][j][k+1]-T[i][j][k])/(z[k+2]-z[k]) ) )

# endif

# define E_CONV  (udTdr_UD1 + vrdTdth_UD1 + wdTdz_UD1)


/* Discretizing the conduction term */
# ifndef  THMCOND_TDEP

#   define d2Tdr2 ( 8*((T[i+1][j][k]-T[i][j][k])/(r[i+2]-r[i])-(T[i][j][k]-T[i-1][j][k])/(r[i+1]-r[i-1]))/(r[i+2]+r[i+1]-r[i]-r[i-1]) )

#   define _rdTdr ( 4*((r[i+1]-r[i-1])*(T[i+1][j][k]-T[i][j][k])/(r[i+2]-r[i]) + (r[i+2]-r[i])*(T[i][j][k]-T[i-1][j][k])/(r[i+1]-r[i-1]))/((r[i+1]+r[i])*(r[i+2]+r[i+1]-r[i]-r[i-1])) )

#   define _r2d2Tdth2 ( 32*((T[i][j+1][k]-T[i][j][k])/(th[j+2]-th[j])-(T[i][j][k]-T[i][j-1][k])/(th[j+1]-th[j-1]) )/((r[i]+r[i+1])*(r[i]+r[i+1])*(th[j+2]+th[j+1]-th[j]-th[j-1])) )

#   define d2Tdz2 ( 8*((T[i][j][k+1]-T[i][j][k])/(z[k+2]-z[k])-(T[i][j][k]-T[i][j][k-1])/(z[k+1]-z[k-1]))/(z[k+2]+z[k+1]-z[k]-z[k-1]) )

#   define  E_COND  ( thmcond0*(d2Tdr2 + _rdTdr + _r2d2Tdth2 + d2Tdz2)/(density0*Cp_WATER) )

# else

#   define COND_TERM_1   ( 4*( thm_cond((T0[i+1][j][k]*(r[i+1]-r[i])+T0[i][j][k]*(r[i+2]-r[i+1]))/(r[i+2]-r[i]))\
         *r[i+1]*(T[i+1][j][k]-T[i][j][k])/(r[i+2]-r[i]) - thm_cond((T0[i][j][k]*(r[i]-r[i-1])+T0[i-1][j][k]*(r[i+1]-r[i]))/(r[i+1]-r[i-1]))\
         *r[i]*(T[i][j][k]-T[i-1][j][k])/(r[i+1]-r[i-1]) )/((r[i+1]+r[i])*(r[i+1]-r[i])) )

#   define COND_TERM_2   ( 8*( thm_cond((T0[i][j+1][k]*(th[j+1]-th[j])+T0[i][j][k]*(th[j+2]-th[j+1]))/(th[j+2]-th[j]))\
         *(T[i][j+1][k]-T[i][j][k])/(th[j+2]-th[j]) - thm_cond((T0[i][j][k]*(th[j]-th[j-1])+T0[i][j-1][k]*(th[j+1]-th[j]))/(th[j+1]-th[j-1]))\
         *(T[i][j][k]-T[i][j-1][k])/(th[j+1]-th[j-1]) )/((r[i]+r[i+1])*(r[i]+r[i+1])*(th[j+1]-th[j])) )

#   define COND_TERM_3   ( 2*( thm_cond((T0[i][j][k+1]*(z[k+1]-z[k])+T0[i][j][k]*(z[k+2]-z[k+1]))/(z[k+2]-z[k]))\
         *(T[i][j][k+1]-T[i][j][k])/(z[k+2]-z[k]) - thm_cond((T0[i][j][k-1]*(z[k+1]-z[k])+T0[i][j][k]*(z[k]-z[k-1]))/(z[k+1]-z[k-1]))\
         *(T[i][j][k]-T[i][j][k-1])/(z[k+1]-z[k-1]) )/(z[k+1]-z[k]) )

#   define  E_COND (COND_TERM_1 + COND_TERM_2 + COND_TERM_3)/(density0*Cp_WATER)


# endif

# define  VHSRC_TERM   ((vhsrc_flag == 'y') ? EXPECTED_T_RISE/(double)EXPECTED_TIME:0)

# define  E_ijk   (VHSRC_TERM + E_COND - E_CONV)


#endif
/* ======================================================================================================================== */
